\documentclass[11pt]{ctexart}
\usepackage{tabularx}
\usepackage{array}
\usepackage{bm}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{algorithm}
\usepackage{graphicx}
\usepackage[table,xcdraw]{xcolor}
%\usepackage{algorithm2e}
\usepackage{algorithmicx}
\usepackage{algpseudocode}
\usepackage{fancyhdr}
\pagestyle{fancy}
\usepackage{color}
\usepackage{xcolor}
\usepackage{float}
\definecolor{keywordcolor}{rgb}{0.8,0.1,0.5}
\usepackage{listings}
\hypersetup{hypertex=true,
	colorlinks=true,
	linkcolor=red,
	anchorcolor=blue,
    citecolor=blue}
\definecolor{mygreen}{rgb}{0,0.6,0}  
\definecolor{mygray}{rgb}{0.5,0.5,0.5}  
\definecolor{mymauve}{rgb}{0.58,0,0.82}  
\lstset{ %  
    backgroundcolor=\color{white},   % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}  
    basicstyle=\footnotesize,        % the size of the fonts that are used for the code  
    breakatwhitespace=false,         % sets if automatic breaks should only happen at whitespace  
    breaklines=true,                 % sets automatic line breaking  
    captionpos=bl,                    % sets the caption-position to bottom  
    commentstyle=\color{mygreen},    % comment style  
    deletekeywords={...},            % if you want to delete keywords from the given language  
    escapeinside={\%*}{*)},          % if you want to add LaTeX within your code  
    extendedchars=true,              % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8  
    frame=single,                    % adds a frame around the code  
    keepspaces=true,                 % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)  
    keywordstyle=\color{blue},       % keyword style  
    %language=Python,                 % the language of the code  
    morekeywords={*,...},            % if you want to add more keywords to the set  
    numbers=left,                    % where to put the line-numbers; possible values are (none, left, right)  
    numbersep=5pt,                   % how far the line-numbers are from the code  
    numberstyle=\tiny\color{mygray}, % the style that is used for the line-numbers  
    rulecolor=\color{black},         % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))  
    showspaces=false,                % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'  
    showstringspaces=false,          % underline spaces within strings only  
    showtabs=false,                  % show tabs within strings adding particular underscores  
    stepnumber=1,                    % the step between two line-numbers. If it's 1, each line will be numbered  
    stringstyle=\color{orange},     % string literal style  
    tabsize=2,                       % sets default tabsize to 2 spaces  
    %title=myPython.py                   % show the filename of files included with \lstinputlisting; also try caption instead of title  
}  
\fancyhf{}
\chead{\textbf{计算方法}}
\fancyhead[r]{\bfseries\thepage}
\fancyhead[l]{\bfseries\rightmark}
\renewcommand{\headrulewidth}{0.4pt} % 注意不用 \setlength
\renewcommand{\footrulewidth}{0pt}

\floatname{algorithm}{算法}
\renewcommand{\algorithmicrequire}{\textbf{输入:}}
\renewcommand{\algorithmicensure}{\textbf{输出:}}
\begin{titlepage}
	\title{\Huge\textbf{计算方法第六、八章作业}\\}
	\author{\Large\textbf{作者}：吴润泽 \and{\Large\textbf{学号}：181860109}\\
		\\
		\and {\Large\textbf{Email}：\href{mailto:181860109@smail.nju.edu.cn}{181860109@smail.nju.edu.cn}}\\}
	\date{\Large\today}
\end{titlepage}
\begin{document}
		\maketitle
	\tableofcontents
	\newpage
	\section*{assignment I}
	\markright{assignment I}
	\addcontentsline{toc}{section}{assignment I}
	\subsection*{problem 1}
	\markright{problem 1}
    \addcontentsline{toc}{subsection}{problem 1}
    \paragraph*{解}设$f(x)=x^2-x-1$。\\
    因为$f(0)=-1<0,f(2)=1>0$，所以$[0,2]$为$f(x)$的有根区间。\\
    根据二分法的误差估计式，误差小于0.05，即$\frac{2}{2^{k+1}}<0.05$，解得$k>4$，
    故至少应二分5次。\\
    因为$f(1)=-1<0$，所以$[1,2]$为$f(x)$的有根区间。\\
    因为$f(1.5)=-0.25<0$，所以$[1.5,2]$为$f(x)$的有根区间。\\
    因为$f(1.75)=\frac{5}{16}>0$，所以$[1.5,1.75]$为$f(x)$的有根区间。\\
    因为$f(1.625)=\frac{1}{64}>0$，所以$[1.5,1.625]$为$f(x)$的有根区间。\\
    因为$f(1.5625)=-\frac{31}{256}>0$，所以$[1.5625,1.625]$为$f(x)$的有根区间。\\
    则取$x^{*}=1.59375$，即为所求。
    \subsection*{problem 3}
	\markright{problem 3}
    \addcontentsline{toc}{subsection}{problem 3}
    \paragraph*{解} 取$x_0$的邻域$[1.3,1.6]$.
    \paragraph*{(1)}
    $\varphi(x)=1+\frac{1}{x^2}$在$[1.3,1.6]$连续,
    $|\varphi^{'}(x)|=|-\frac{2}{x^3}|$，\\
    $|\varphi^{'}(1.5)|=\frac{2}{1.5^3}\approx0.593<1$，
    故迭代$x_{k+1}=1+\frac{1}{x^2_k}$在$[1.3,1.6]$局部收敛。
    \paragraph*{(2)}
    $\varphi(x)=(1+x^2)^{1/3}$在$[1.3,1.6]$连续,
    $|\varphi^{'}(x)|=\frac{2}{3}|\frac{x}{(1+x^2)^{2/3}}|$，\\
    $|\varphi^{'}(1.5)|=\sqrt[3]{\frac{16}{169}}\approx0.456<1$，
    故迭代$x_{k+1}=(1+x_k^2)^{1/3}$在$[1.3,1.6]$局部收敛。
    \paragraph*{(3)}
    $\varphi(x)=\frac{1}{\sqrt{x-1}}$在$[1.3,1.6]$连续,
    $|\varphi^{'}(x)|=\frac{1}{2}|(x-1)^{-\frac{3}{2}}|$，\\
    $|\varphi^{'}(1.5)|=\sqrt{2}\approx1.414>1$，
    故迭代$x_{k+1}=\frac{1}{\sqrt{x_k-1}}$发散。

    由于(2)中的$|\varphi^{'}(1.5)|$较小，故取(2)得迭代公式计算，
    其中$$|\varphi^{'}(x)|<\frac{2}{3}\frac{1.6}{(1+1.3^2)^{2/3}}\approx0.522=L<1,$$
    若使结果具有四位有效数字，则
    $$|x_k-x^{*}|\leq\frac{L}{1-L}|x_k-x_{k-1}|<\frac{1}{2}\times10^{-3},$$
    可推得
    $$|x_k-x_{k-1}|<\frac{1-L}{L}\times\frac{1}{2}\times10^{-3}<\frac{1}{2}\times10^{-3}.$$
    \begin{tabular}{|c|c|c|c|c|c|}
        \hline
        \textbf{k} & \textbf{$x_k$} & \textbf{k} & \textbf{$x_k$} & \textbf{k} & \textbf{$x_k$} \\ \hline
        1 & 1.481\ 248 & 3 & 1.468\ 817 & 5 & 1.466\ 243 \\ \hline
        2 & 1.472\ 706 & 4 & 1.467\ 048 & 6 & 1.465\ 877 \\ \hline
    \end{tabular}\\
    由于$|x_6-x_5|<\frac{1}{2}\times10^{-3}$，故可取$x^{*}\approx x_6=1.466$.
    \subsection*{problem 7}
	\markright{problem 7}
    \addcontentsline{toc}{subsection}{problem 7}
    \paragraph*{(1)}牛顿法
    $$\forall x\in[1,2],f(1)<0,f(2)>0,f^{'}(x)=3(x^2-1)\geq0,f^{''}(x)=6x>0.$$
    $$x_{k+1}=x_k-\frac{x_k^2-3x_k-1}{3x_k^2-3}=\frac{2x^2_k+1}{3(x^2_k-1)},k=0,1,2,\cdots$$
    计算得$$x_1=1.888\ 888\ 89\ ,\ x_2=1.879\ 451\ 67,$$
    因为$|x_2-x^{*}|<\frac{1}{2}\times10^{-3}$，所以$x^{*}\approx x_2=1.879.$
    \paragraph*{(2)}弦截法
    $$x_{k+1}=x_k-\frac{(x_k-x_{k-1})f(x_k)}{f(x_k)-f(x_{k-1})}=\frac{x^2_k x_{k-1}+x_k x_{k-1}^2+1}{x^2_k+x_k x_{k-1}+x^2_{k-1}-3},k=1,2,\cdots$$
    计算得
    $$x_2=1.881\ 093\ 94,\ x_3=1.879\ 411\ 06,$$
    因为$|x_3-x^{*}|<\frac{1}{2}\times10^{-3}$，所以$x^{*}\approx x_3=1.879.$
    \paragraph*{(3)}抛物线法
    $$\left\{
        \begin{aligned}
            &x_{k+1}=x_k-\frac{2f(x_k)}{w\pm\sqrt{w^2-4f(x_k)f[x_k,x_{k-1},x_{k-2}]}},\\
            &w=f[x_k,x_{k-1}]+f[x_k,x_{k-1},x_{k-2}](x_k-x_{k-1}),
        \end{aligned}\right .
    $$
    $\begin{aligned}
    &f(x_0)=-3,f(x_1)=17,f(x_2)=1,f[x_0,x_1]=10,f[x_2,x_1]=16,\\
    &f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}=6,w=16+6(2-3)=10,\\
\end{aligned}\\$
    计算得$$x_3=1.893\ 149\ 82,x_4=1.879\ 135\ 26,$$
    因为$|x_4-x^{*}|<\frac{1}{2}\times10^{-3}$，所以$x^{*}\approx x_4=1.879.$
    \subsection*{problem 12}
	\markright{problem 12}
    \addcontentsline{toc}{subsection}{problem 12}
    \paragraph*{解} $f(x)=x^3-a$，故$f^{'}(x)=3x^2$，牛顿法迭代公式为
    $$x_{k+1}=x_k-\frac{f(x_k)}{f^{'}(x_k)}=x_k-\frac{3x_k^3-a}{3x_k^2}=\frac{2x^3_k+a}{3x^2_k}=\varphi(x_k),k=0,1,2,\cdots.$$
    当$a\neq0$时，$\sqrt[3]{a}$为$f(x)=0$的单根，$\varphi^{'}(x^{*})=\frac{2({x^{*}}^3-a)}{3{x^{*}}^3}=0$，$\varphi^{''}(x^{*})=\frac{2a}{{x^{*}}^4}\neq 0$，
    迭代公式为$x_{k+1}=\frac{2x_k^3+a}{3x_k^2}$，在$x^{*}$附近是平方收敛的。

    \noindent 当$a=0$时，迭代公式退化为$x_{k+1}=\frac{2}{3}x_k=\varphi(x_k)$，$\varphi(x^{*})=0,\varphi^{'}(x^{*})\neq 0$，
    线性收敛。
    \subsection*{problem 13}
	\markright{problem 13}
    \addcontentsline{toc}{subsection}{problem 13}
    \paragraph*{解} $x^{*}=\sqrt{a},a\neq 0$为方程$f(x)$的单根，$f^{'}(x)=\frac{2a}{x^3}$，牛顿法迭代公式为
    $$x_{k+1}=x_k-\frac{f(x_k)}{f^{'}(x_k)}=x_k-\frac{x_k^3-ax_k}{2a}=\frac{1}{2a}(3ax_k-x_k^3)=\varphi(x_k),k=0,1,2,\cdots.$$
    $\varphi^{'}(x^{*})=\frac{1}{2a}(3a-3a{x^{*}}^2)=0,\varphi^{''}(x^{*})=3x^{*}\neq 0$。

\noindent 因此迭代公式为$x_{k+1}=\frac{1}{2a}(3ax_k-x_k^3)$，在$x^{*}$附近是平方收敛的。
    \newpage
    \section*{assignment II}
	\markright{assignment II}
	\addcontentsline{toc}{section}{assignment II}
	\subsection*{problem 1}
	\markright{problem 1}
    \addcontentsline{toc}{subsection}{problem 1}
    \paragraph*{(1)}系数矩阵为
    $$
    \begin{aligned}
        A=
        \begin{bmatrix}
        5&2&1\\-1&4&2\\2&-3&10\\    
        \end{bmatrix}
    \end{aligned}
    $$较易验证其满足$$a_{ii}>\sum \limits_{j=1,j\neq i}^{n}|a_{ij}|,\ (i=1,2,\cdots,n).$$
    即系数矩阵严格对角占优，故Jacobi迭代法和Gauss-Seidel迭代法均收敛。
    \paragraph*{(2)}\textbf{Jacobi迭代法}
    $$\left\{
        \begin{aligned}
            &x_1^{(k+1)}=-\frac{2}{5}x_2^{(k)}-\frac{1}{5}x_3^{k}-\frac{12}{5},\\
            &x_2^{(k+1)}=\frac{1}{4}x_1^{(k)}-\frac{1}{2}x_3^{(k)}+5,\\
            &x_3^{(k+1)}=-\frac{1}{5}x_1^{(k)}+\frac{3}{10}x_2^{(k)}+\frac{3}{10}.\\
        \end{aligned}
        \right.
    $$
    取$x^{(0)}=(1,1,1)^{T}$，则迭代17次可达到精度要求，即
    $$x^{(17)}=(-4.000\ 03,\ 2.999\ 98,\ 2.000\ 03)^{T}.$$
    \textbf{Gauss-Seidel迭代法}
    $$\left\{
        \begin{aligned}
            &x_1^{(k+1)}=-\frac{2}{5}x_2^{(k)}-\frac{1}{5}x_3^{k}-\frac{12}{5},\\
            &x_2^{(k+1)}=\frac{1}{4}x_1^{(k+1)}-\frac{1}{2}x_3^{(k)}+5,\\
            &x_3^{(k+1)}=-\frac{1}{5}x_1^{(k+1)}+\frac{3}{10}x_2^{(k+2)}+\frac{3}{10}.\\
        \end{aligned}
        \right.
    $$
    取$x^{(0)}=(1,1,1)^{T}$，则迭代8次可达到精度要求，即
    $$x^{(8)}=(-4.000\ 02,\ 2.999\ 99,\ 2.000\ 00)^{T}.$$
    \subsection*{problem 5}
	\markright{problem 5}
    \addcontentsline{toc}{subsection}{problem 5}
    \paragraph*{(1)}
    $
    \begin{aligned}
        &D=
        \begin{bmatrix}
        1&&\\&1&\\&&1\\    
        \end{bmatrix},
        L= \begin{bmatrix}
            0&&\\-0.4&0&\\-0.4&-0.8&0\\    
            \end{bmatrix},
        &U= \begin{bmatrix}
            0&-0.4&-0.4\\&0&-0.8\\&&0\\    
            \end{bmatrix},
    \end{aligned}$\\
    \textbf{Jacobi}迭代矩阵为
    $$
    \begin{aligned}
        &B_{Jacobi}=D^{-1}(L+U)=
        \begin{bmatrix}
        0&-0.4&-0.4\\-0.4&0&-0.8\\-0.4&-0.8&0\\    
        \end{bmatrix},\\
    \end{aligned}
    $$
    $\lambda=(-1.09,0.29,0.80)\rightarrow \rho(B_{Jacobi})=1.09>1$
    所以Jacobi迭代法不收敛。
    \textbf{Gauss-Seidel}迭代矩阵为
    $$
    \begin{aligned}
        &B_{G-S}=(D-L)^{-1}U=
        \begin{bmatrix}
        0&-0.4&-0.4\\0&0.16&-0.64\\0&0.032&0.672\\    
        \end{bmatrix},
    \end{aligned}$$
    $\lambda=(0,0.20,0.63)\rightarrow \rho(B_{G-S})=0.63<1$
    所以Gauss-Seidel迭代法收敛。
    \paragraph*{(2)}
    $
    \begin{aligned}
        &D=
        \begin{bmatrix}
        1&&\\&1&\\&&1\\    
        \end{bmatrix},
        L= \begin{bmatrix}
            0&&\\-1&0&\\-2&-2&0\\    
            \end{bmatrix},
        &U= \begin{bmatrix}
            0&-2&2\\&0&-1\\&&0\\    
            \end{bmatrix},
    \end{aligned}$\\
    \textbf{Jacobi}迭代矩阵为
    $$
    \begin{aligned}
        &B_{Jacobi}=D^{-1}(L+U)=
        \begin{bmatrix}
        0&-2&2\\-1&0&-1\\-2&-2&0\\    
        \end{bmatrix},
    \end{aligned}$$
    $\lambda=(0,0,0)\rightarrow \rho(B_{Jacobi})=0<1$
    所以Jacobi迭代法收敛。\\
    \textbf{Gauss-Seidel}迭代矩阵为
    $$
    \begin{aligned}
        &B_{G-S}=(D-L)^{-1}U=
        \begin{bmatrix}
        0&-2&-2\\0&2&-3\\0&0&2\\    
        \end{bmatrix},
    \end{aligned}$$
    $\lambda=(0,2,2)\rightarrow \rho(B_{G-S})=2>1$
    所以Gauss-Seidel迭代法不收敛。
    \subsection*{problem 9}
	\markright{problem 9}
    \addcontentsline{toc}{subsection}{problem 9}
    SOR迭代法公式为
    $$\left\{
        \begin{aligned}
            &x_1^{(k+1)}=x_1^{(k)}+\omega(\frac{1}{4}-x_1^{(k)}+\frac{1}{4}x_2^{(k)}),\\
            &x_2^{(k+1)}=x_2^{(k)}+\omega(1+\frac{1}{4}x_1^{(k+1)}-x_2^{(k)}+\frac{1}{4}x_3^{(k)}),\\
            &x_3^{(k+1)}=x_3^{(k)}+\omega(-\frac{3}{4}+\frac{1}{4}x_2^{(k+1)}-x_3^{(k)}).\\
        \end{aligned}
        \right.
    $$
    当取$\omega=1.03$，初值$x^{(0)}=(0,0,0)^{T}$时，迭代5次可达到精度要求，即
    $$x^{(5)}=(0.500\ 004\ 5,\ 1.000\ 001\ 6,\ -4.999\ 997)^{T}.$$
    当取$\omega=1$，初值$x^{(0)}=(0,0,0)^{T}$时，迭代6次可达到精度要求，即
    $$x^{(5)}=(0.500\ 003\ 8,\ 1.000\ 001\ 9,\ -4.999\ 995)^{T}.$$
    当取$\omega=1.1$，初值$x^{(0)}=(0,0,0)^{T}$时，迭代6次可达到精度要求，即
    $$x^{(5)}=(0.500\ 003\ 6,\ 0.999\ 998\ 5,\ -5.000\ 000\ 0)^{T}.$$
    \subsection*{problem 14}
	\markright{problem 14}
    \addcontentsline{toc}{subsection}{problem 14}
    \paragraph*{证明} 
    由A的顺序主子式
    $\triangle_2=\begin{bmatrix}
        1&a\\a&1\\
    \end{bmatrix}
    =1-a^2>0$，得$-1<a<1$，
    而$\triangle_3=det(A)=2a^3-2a^2+1=(a-1)^2(2a+1)>0$，得$a>-\frac{1}{2}$，
    即$-\frac{1}{2}<a<1$时，A正定。\\
    对于Jacobi迭代矩阵为
    $$J=\begin{bmatrix}
        0&-a&-a\\-a&0&-a\\-a&-a&0\\
    \end{bmatrix}$$
    $det(\lambda I-J)=(\lambda-a)^2(\lambda+2a)=0$，$\rho(J)=|2a|<1$，
    即$-\frac{1}{2}<a<\frac{1}{2}$时Jacboi迭代法收敛。综上得证。
\newpage
\section*{附录}
\markright{附录}
\addcontentsline{toc}{section}{附录}
\subsection*{牛顿法}
\markright{牛顿法}
\addcontentsline{toc}{subsection}{牛顿法}
\begin{lstlisting}[language=matlab,title=newton.m]
function [x] = newton(fun,x0,x_,e)
    %fun为目标函数
    %x0为牛顿法起点
    %x_为精确解，e为x-x_的允许误差
    x=x0;k=0;
    syms a;y=fun(a);
    dfun=matlabFunction(diff(y));
    while abs(x-x_)>e
        x0=x;
        x=x0-feval(fun,x0)/feval(dfun,x0);
        k=k+1;X(k)=x;K(k)=k;
    end
    A=[K;X];fprintf("%d %.8f\n",A);
end
\end{lstlisting}
\subsection*{弦截法}
\markright{弦截法}
\addcontentsline{toc}{subsection}{弦截法}
\begin{lstlisting}[language=matlab,title=secant.m]
function [x] = secant(fun,x0,x1,x_,e)
    %fun为目标函数
    %x0,x1为弦截法起点
    %x_为精确解，e为x-x_的允许误差
    x=x1;y=x0;k=0;
    while abs(x-x_)>e
        z=x-(feval(fun,x)*(x-y))/(feval(fun,x)-feval(fun,y));
        y=x;x=z;
        k=k+1;X(k)=x;K(k)=k;
    end
    A=[K;X];fprintf("%d %.8f\n",A);
end
\end{lstlisting}
\newpage
\subsection*{抛物线法}
\markright{抛物线法}
\addcontentsline{toc}{subsection}{抛物线法}
\begin{lstlisting}[language=matlab,title=parobola.m]
function [x] = parobola(fun,x0,x1,x2,x_,e)
    %fun为目标函数
    %x0,x1,x2为抛物线法选择的三个点
    %x_为精确解，e为x-x_的允许误差
    k=0;
    z=x0;
    y=x1;
    x=x2;
    while abs(x-x_)>e
        h1=y-z;
        h2=x-y;
        c1=(feval(fun,y)-feval(fun,z))/h1;
        c2=(feval(fun,x)-feval(fun,y))/h2;
        d=(c2-c1)/(h2+h1);
        w=c2+h2*d;
        xi=x-(2*feval(fun,x))/(w+sign(w)*\
           sqrt(w^2-4*feval(fun,x)*d));
        z=y;
        y=x;
        x=xi;
        k=k+1;X(k)=x;K(k)=k;
    end
    A=[K;X];fprintf("%d %.8f\n",A);
end
\end{lstlisting}
\newpage
\subsection*{Jacobi迭代法}
\markright{Jacobi迭代法}
\addcontentsline{toc}{subsection}{Jacobi迭代法}
\begin{lstlisting}[language=matlab,title=LinJacobi.m]
function [x,k] = LinJacobi(A,b,x0,tol)
    %A为系数矩阵，b为矩阵右侧
    %x0为初始x，t0为x(k+1)-x(k)的精度
    %x为方程组的解，k为迭代次数
    max=3000;
    D=diag(diag(A));
    L=-tril(A,-1);
    U=-triu(A,1);
    B=D\(L+U);
    f=D\b;
    x=B*x0+f;
    k=1;
    while norm(x-x0,Inf)>=tol
        x0=x;
        x=B*x0+f;
        k=k+1;
        if(k>=max)
            disp('超出迭代次数');
            break;
        end
    end
end
\end{lstlisting}
\newpage
\subsection*{Gauss-Seidel迭代法}
\markright{Gauss-Seidel迭代法}
\addcontentsline{toc}{subsection}{Gauss-Seidel迭代法}
\begin{lstlisting}[language=matlab,title=LinGauSeid.m]
function [x,k] = LinGauSeid(A,b,x0,tol)
    %A为系数矩阵，b为矩阵右侧
    %x0为初始x，t0为x(k+1)-x(k)的精度
    %x为方程组的解，k为迭代次数
    max=3000;
    D=diag(diag(A));
    L=-tril(A,-1);
    U=-triu(A,1);
    G=(D-L)\U;
    f=(D-L)\b;
    x=G*x0+f;
    k=1;
    while norm(x-x0,Inf)>=tol
        x0=x;
        x=G*x0+f;
        k=k+1;
        if(k>=max)
            disp('超出迭代次数');
            break;
        end
    end
end
\end{lstlisting}
\newpage
\subsection*{SOR迭代法}
\markright{SOR迭代法}
\addcontentsline{toc}{subsection}{SOR迭代法}
\begin{lstlisting}[language=matlab,title==LinSOR.m]
function [x,k] = LinSOR(A,b,w,x0,cut,x_)
    %A为系数矩阵，b为矩阵右侧，x_为精确解
    %x0为初始x，cut为x-x_的精度
    %x为方程组的解，k为迭代次数
    max=30000;
    D=diag(diag(A));
    L=-tril(A,-1);
    U=-triu(A,1);
    G=(D-w.*L)\((1-w).*D+w.*U);
    f=(D-w.*L)\(w.*b);
    x=G*x0+f;
    k=1;
    while norm(x-x_)>=cut
        x0=x;
        x=G*x0+f;
        k=k+1;
        if(k>=max)
            disp('超出迭代次数');
            break;
        end
    end
end\end{lstlisting}
\end{document}